library(DBI)
library(sf)
library(dplyr)
library(gimms)
library(raster)
library(ncdf4)
library(rgdal)
library(spdep)

library(stringr)
library(ggplot2)
library(cowplot)
library(spdep)
library(tibble)
library(RColorBrewer)

library(kableExtra)

library(gridExtra)
library(grid)
library(png)

source("f_plotspatial.R")
memory.limit(size = 160000)
## [1] 160000
compute.effect.size <- function(df.coef){

df.es <- df.coef[,1:3]
n <- length(df.es$Parameter)
m <- df.es$Parameter


for(i in 1:n){
  if(grepl("log",m[i])==T){
    df.es[i,2] <- 1.01^(df.coef[i,2])
    df.es[i,3] <- 1.01^(df.coef[i,2])*log(1.01)*df.coef[i,3]
  }else if(m[i] %in% c("NDVI","Bog","Arable","Forest","Temp")){
    df.es[i,2] <- exp(df.coef[i,2]*0.01)
    df.es[i,3] <- 0.01*df.coef[i,2]*df.coef[i,3]
  }else if(m[i] %in% c("TNdep")){
    df.es[i,2] <- exp(df.coef[i,2]*25)
    df.es[i,3] <- 25*exp(df.coef[i,2]*25)*df.coef[i,3]
  }else if(m[i] %in% c("TSdep")){
    df.es[i,2] <- exp(df.coef[i,2]*11)
    df.es[i,3] <- 11*exp(df.coef[i,2]*11)*df.coef[i,3]
  }else {
    df.es[i,2] <- exp(df.coef) # effect size for temperature???
  }
}

df.es$Percent <- 100*(df.es$Estimate-1)

return(df.es)

}

# other es
eslog <- 1.01
es <- 0.01
library(captioner)
fig_nums <- captioner() 
table_nums <- captioner(prefix = "Table")

1 Catchment data for Norway, Sweden and Finland

1.1 Catchment polygons

The catchment and TOC data are stored on the PostgreSQL database ecco_biwa, accessed through Rstudio using the st_read function (sf package Pebesma (2021)) and the dplyr package (Wickhman2021?). Each catchment polygon is linked to the corresponding TOC concentration via an “ebint” (identifyer), internal to the ecco_biwa database.

The catchment polygons were designed by an elevation model and are assigned to the studied lake based on the distance to the sampling coordinates.

con <- dbConnect(RPostgreSQL::PostgreSQL(),user = "camille.crapart", password = "camille",host = "vm-srv-wallace.vm.ntnu.no", dbname = "nofa")
#catchment.geom <- tbl(con,sql("SELECT ebint, geom FROM catchments.lake_catchments"))

ebint.tbl <- tbl(con,sql("SELECT ebint FROM catchments.lake_catchments")) 
ebint.catch <- pull(ebint.tbl,ebint)

ebint.waterchem <- tbl(con,sql("SELECT ebint,nation FROM environmental.north_euro_lake_surv_1995")) %>% tbl_df()

ebint <- intersect(ebint.waterchem$ebint, ebint.catch)

country <- filter(ebint.waterchem, ebint %in% ebint)
saveRDS(country,"country.Rdata")
catchment.poly.1000 <- st_read(con,query = "SELECT ebint, geom FROM catchments.lake_catchments WHERE ebint IN (SELECT ebint FROM environmental.north_euro_lake_surv_1995)")
saveRDS(catchment.poly.1000,"catchment.poly.1000.Rdata")

The catchment polygons were reprojected in WGS84 to match the satellite data that will be extracted later, and in EPSG:3035 (European projection) to match the CORINE land cover data.

catchment.poly <- st_transform(catchment.poly.1000, CRS("EPSG:4326"))
saveRDS(catchment.poly,"catchment.poly.Rdata")
catchment.poly.corine <- st_transform(catchment.poly.1000, CRS("EPSG:3035"))
saveRDS(catchment.poly.corine, "Fennoscandia_NTNU/catchment.poly.Rdata")

1.2 TOC

The TOC data was measured after the Northern Lake Survey 1995 (Henriksen et al. 1998). It was stored in the ecco_biwa database (https://github.com/andersfi/BiWa_ECCO).

toc.tbl <- tbl(con,sql("SELECT ebint,toc_mg_l,longitude,latitude,dist_closest_ebint,dist_2nd_closest_ebint FROM environmental.north_euro_lake_surv_1995"))
toc.df <- as.data.frame(toc.tbl)
saveRDS(toc.df,"toc.df.Rdata")

1.3 NDVI

NDVI values are extracted from the GIMMS NDVI3g dataset (TheNationalCenterForAtmosphericResearch2018?), stored on http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, in the “ecocast” dataset and accessed via the “gimms” package (Destch2021?) on Rstudio.

Data is taken bi-monthly (Tucker et al. 2005). Raster layers of all slices are downloaded using the gimmsdownload function, then monthly composite are calculated using the monthlyComposites function. The maximum NDVI value is taken for each pixel. Afterwards the NDVI values for each catchment polygon are extracted using the extract function from the raster package (Hijmans et al. 2018).

The values stored in the ecocast dataset are composed of the NDVI value and a flag value indicating the goodness of the data (Detsch 2021). All the NDVI values extracted, corresponding to the values for the studied catchments, had a flag of 1, indicating a good value. The NDVI value was retrieved from the stored value using the formula \(floor(ndvi3g/10)/1000\) (Detsch 2021), and then the mean of the 3 summer values (June, July and August) was computed for each catchment.


# Download Gimms NDVI

ndvi.1994 <- downloadGimms(x= 1994,y=1994,dsn = "NDVI")
ndvi.max <- monthlyComposite(ndvi.1994,monthlyIndices(ndvi.1994))
saveRDS(ndvi.max,"ndvi.max.rds")

# Load data

ndvi.max <- readRDS("ndvi.max.rds")
catchment.poly <- readRDS("catchment.poly.Rdata")

# Computes summer mean
summer.mean <- raster::stack(ndvi.max[[6]],ndvi.max[[7]],ndvi.max[[8]]) %>% mean()

# Crops raster to Scandinavia
summer.scandinavia <- raster::crop(summer.mean,c(0,35,55,73))

# Re-introduces NA
summer.scan <- reclassify(summer.scandinavia, cbind(-Inf, 0, NA), right=FALSE)
saveRDS(summer.scan, "summer.scan.rds")

# Extract NDVI for each catchment
summer.ndvi <- raster::extract(summer.mean,catchment.poly, fun = mean, df = T, sp = T)
names(summer.ndvi) <- c("ebint","ndvi")
summer.ndvi$ndvi.value <- floor(summer.ndvi$ndvi/10)/1000
summer.ndvi$flag.value <- summer.ndvi$ndvi - floor(summer.ndvi$ndvi/10)*10 + 1 

# Saves NDVI file
saveRDS(summer.ndvi, "summer.ndvi.Rdata")
write.csv(summer.ndvi,"summer.ndvi.csv")
summer.scan <- readRDS("summer.scan.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")

ndvi.plot <- (floor(summer.scan)/10)/1000

ndvi.plot.spdf <- ndvi.plot %>% as("SpatialPixelsDataFrame") %>% spTransform(projection(wr.sf.95)) %>% as.data.frame()

map.ndvi <- ggplot() + geom_tile(data = ndvi.plot.spdf, aes(x = x, y = y, fill = layer, width = 10000, height = 10000))+
  scale_fill_distiller(type = "seq", palette = 2, direction = 1)+labs(fill = "Summer NDVI 1995", subtitle = "a)")+
  theme_void(base_size = 8)+
  theme(legend.position = "bottom")
ggsave(plot = map.ndvi, filename = "NDVI/map.ndvi.nsf.png", dpi = "retina", units = "cm")
knitr::include_graphics("NDVI/map.ndvi.nsf.png")

fig_nums(map_ndvi, “Map of summer NVDI in Norway, Sweden and Finland in 1994”)

1.4 Runoff

The runoff data was downloaded from the CORDEX database CORDEX (2021). The data is available at https://esg-dn1.nsc.liu.se/projects/esgf-liu/.

The names of the variables in the CORDEX project follow the CF Metadata convention: http://cfconventions.org/

THe mean of the 30 years between 1970 and 2000 is used (p11 on EURO-CORDEX guidelines https://www.euro-cordex.net/imperia/md/content/csc/cordex/guidance_for_euro-cordex_climate_projections_data_use__2021-02_1_.pdf). This time period matches the temperature and precipitation data provided by WorldClim.

runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)

runoff.stack <- raster::stack(runoff.raster, bands = c(1441:1812)) # Correspond to January 1970 to December 2000
runoff.mean <- runoff.stack %>% mean() 

catchment.poly <- readRDS("catchment.poly.Rdata")
runoff.fennoscandia <- raster::extract(runoff.mean, catchment.poly, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.fennoscandia) <- c("ebint","runoff")
saveRDS(runoff.fennoscandia,"CORDEX/runoff.fennoscandia.rds")

1.5 Temperature, precipitation and slope

Temperature and precipitation are downloaded from the database Worldclim. Temperature is the mean annual temperature from 1970 to 2000 in C, or bio1 in the bioclimatic dataset. The original data is in C*10 so the data is divided by 10 before use. Precipitation is the mean annual precipitation in mm or bio12 in the bioclimatic dataset : https://www.worldclim.org/data/worldclim21.html

bio.fennoscandia <- raster::getData(name = "worldclim", var = "bio", res = 2.5)

annual.temp <- raster::extract(bio.fennoscandia[[1]], catchment.poly, fun = mean, df = T, sp = T)
saveRDS(annual.temp,"annual.temp.Rdata")

annual.prec <- raster::extract(bio.fennoscandia$bio12, catchment.poly, fun = mean, df = T, sp = T)
saveRDS(annual.prec,"annual.prec.Rdata")

The altitude data was accessed through the “getData” function, and comes from the NASA digital elevation model (SRTM 90 m). The slope was then computed using the “terrain” function from the raster package.

alt.nor <- getData("alt",country = "NOR", mask = T)
alt.swe <- getData("alt",country = "SWE", mask = T)
alt.fin <- getData("alt",country = "FIN", mask = T)

alt.list <- list(alt.nor,alt.swe,alt.fin)
alt <- do.call(raster::merge,alt.list)

alt.fennoscandia <- extract(alt,remote.set[,c("longitude","latitude")])
alt.df <- data.frame(ebint = remote.set$ebint, alt = alt.fennoscandia )
saveRDS(alt.df, "alt.df.Rdata")

slope.nor <- terrain(alt.nor, unit = "degrees")
slope.swe <- terrain(alt.swe, unit = "degrees")
slope.fin <- terrain(alt.fin, unit = "degrees")
slope.list <- list(slope.nor,slope.swe,slope.fin)
slope <- do.call(raster::merge,slope.list)

slope.nona <- reclassify(slope,cbind(NA,100), right = NA) #reintroduces NA

slope.fennoscandia <- raster::extract(slope,catchment.poly, sp = T, df = T, fun = mean)
saveRDS(slope.df, "slope.df.Rdata")

1.6 Land cover

The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/lcc-2000-2006. The previous version (1995-2000) excludes Norway so the 2000 version was preferred.

knitr::include_graphics("Corine_Land_Cover/Map_clc.png")

fig_nums(map_clc, “Map of Land Cover in Norway, Sweden and Finland in 2000”)

The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:

Arable land: * 12: non-irrigated arable land * 16: fruit trees and berries plantations * 18: pastures * 19: annual crops associated with permanent crops * 20: Complex cultivation patterns

Bogs: * 36: peat bogs * 35: inland marshes, excluded because all zeros in the studied area.

Forest: * 23: Broad-leaved forest * 24: Coniferous forest * 25: Mixed forest

Bare: 31: Bare rocks 32: Sparsely vegetated areas 33: Burnt areas

corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
clc.remote <- raster::extract(corine,catchment.poly.corine)

saveRDS(clc.remote,"clc.remote.Rdata")

The proportion of each land-cover catergory was computed for each catchment.

clc.remote <- readRDS("clc.remote.Rdata")
clc.tab <- sapply(clc.remote, function(x) tabulate(x,45))
clc.tab.area <- clc.tab*prod(res(corine))
catchment.area <- colSums(clc.tab.area)
clc.tab.prop <- sweep(clc.tab.area,2,catchment.area, FUN = "/")
names(clc.tab.prop) <- catchment.poly.corine$ebint
saveRDS(clc.tab.prop,"clc.tab.prop.Rdata")

Then, a dataframe for the category of interest (bogs, arable land, forest, bare land) was created and saved.

clc.tab.prop <- readRDS("clc.tab.prop.Rdata") %>% as.data.frame()

bogs <- clc.tab.prop[36,] %>% t() %>% as.data.frame() %>% setNames("bogs") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(bogs,"bogs.rds")
arable <- colSums(clc.tab.prop[c(12,16,18,19,20),]) %>% as.data.frame() %>% setNames("arable") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(arable,"arable.rds")
forest <- colSums(clc.tab.prop[c(23,24,25),]) %>% as.data.frame() %>% setNames("forest") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(forest,"forest.rds")
bare <- colSums(clc.tab.prop[c(31,32,33),]) %>% as.data.frame() %>% setNames("bare") %>% tibble::rownames_to_column(var = "ebint")
saveRDS(bare,"bare.rds")

1.7 Nitrogen and sulfur deposition

The nitrogen and sulfur deposition data were extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2000 was used as it is the first model available.

N-deposition (NOx and NH3) is the sum of: * Dry deposition of oxidized nitrogen per m2 grid DDEP_OXN_m2Grid, in mg/m2 * Wet deposition of oxidized nitrogen WDEP_OXN, in mg/m2 * Dry deposition of oxidized nitrogen per m2 grid DDEP_RDN_m2Grid, in mg/m2 * Wet deposition of reduced nitrogen WDEP_RDN, in mg/m2

EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
catchment_poly <- readRDS("catchment.poly.Rdata")

woxn <- raster(EMEP_file, varname = "WDEP_OXN") 
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")

woxn_df <- extract(woxn,catchment.poly, sp = T, df = T, na.rm = T, fun = mean) 
names(woxn_df) <- c("ebint","woxn")
doxn_df <- extract(doxn,catchment.poly, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_df) <- c("ebint","doxn")
wrdn_df <- extract(wrdn,catchment.poly, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_df) <- c("ebint","wrdn")
drdn_df <- extract(drdn,catchment.poly, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_df) <- c("ebint","drdn")

ndep.df <- merge(woxn_df,doxn_df, by = "ebint") %>% merge(wrdn_df, by = "ebint") %>% merge(drdn_df, by = "ebint")
ndep.df$tndep <- rowSums(ndep.df@data[,c(2:5)])

saveRDS(ndep.df,"ndep.df.Rdata")

The S-deposition is composed of the sum of: * Dry deposition of oxidized sulphur per m2 grid DDEP_SOX_m2Grid, in mg/m2 * Wet deposition of oxidized sulphur WDEP_SOX, in mg/m2

library(ncdf4)
EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
catchment_poly <- readRDS("catchment.poly.Rdata")

wsox <- raster(EMEP_file, varname = "WDEP_SOX") 
dsox <- raster(EMEP_file, varname = "DDEP_SOX_m2Grid")

wsox_df <- extract(wsox,catchment_poly, sp = T, df = T, na.rm = T, fun = mean) 
names(wsox_df) <- c("ebint","wsox")

dsox_df <- extract(dsox,catchment_poly, sp = T, df = T, na.rm = T, fun = mean)
names(dsox_df) <- c("ebint","dsox")

sdep.df <- merge(wsox_df,dsox_df, by = "ebint")
sdep.df$tsdep <- rowSums(sdep.df@data[,c(2:3)])

saveRDS(sdep.df,"sdep.df.rds")

1.8 Gathering and cleaning the data

The data was gathered and merged into a single dataframe, based on the catchment identifyer (“ebint”).

country <- readRDS("country.Rdata")
toc.df <- readRDS("toc.df.Rdata")
summer.ndvi <- readRDS("summer.ndvi.Rdata")
runoff.fennoscandia <- readRDS("CORDEX/runoff.fennoscandia.rds")

annual.temp <- readRDS("annual.temp.Rdata")
names(annual.temp) <- c("ebint","temp")

annual.prec <- readRDS("annual.prec.Rdata")
names(annual.prec) <- c("ebint","prec")

slope.fennoscandia <- readRDS("slope.fennoscandia.Rdata")
names(slope.fennoscandia) <- c("ebint","slope")

alt.df <- readRDS("alt.df.Rdata")

bogs <- readRDS("bogs.rds")
arable <- readRDS("arable.rds")
forest <- readRDS("forest.rds")
bare <- readRDS("bare.rds")

ndep <- readRDS("ndep.df.Rdata")
sdep <- readRDS("sdep.df.rds")

fennoscandia_raw <- merge(country,toc.df, by = "ebint") %>% 
  merge(summer.ndvi, by = "ebint") %>% 
  merge(runoff.fennoscandia,by = "ebint") %>% 
  merge(annual.temp, by ="ebint") %>% 
  merge(annual.prec, by = "ebint") %>%
  merge(slope.fennoscandia, by = "ebint") %>%
  merge(alt.df, by = "ebint") %>%
  merge(bogs, by = "ebint") %>%
  merge(arable, by = "ebint")%>%
  merge(forest,by = "ebint") %>%
  merge(bare, by = "ebint") %>%
  merge(ndep, by = "ebint") %>% 
  merge(sdep, by = "ebint")

fennoscandia_raw$uncertain <- ifelse(fennoscandia_raw$dist_closest_ebint/fennoscandia_raw$dist_2nd_closest_ebint > 0.5, FALSE, TRUE)
fennoscandia_raw$uncertain[which(is.na(fennoscandia_raw$uncertain)==TRUE)] <- FALSE

fennoscandia_raw <- fennoscandia_raw %>% filter(toc_mg_l >= 0)
fennoscandia_raw <- fennoscandia_raw %>% filter(runoff > 0)
fennoscandia_raw <- fennoscandia_raw %>% filter(ndvi.value > 0)
fennoscandia_raw <- fennoscandia_raw %>% filter(is.na(slope) == F)
fennoscandia_raw <- fennoscandia_raw %>% filter(alt != 0)
fennoscandia_raw <- fennoscandia_raw %>% filter(uncertain == FALSE)
fennoscandia_raw$uncertain <- NULL

saveRDS(fennoscandia_raw,"fennoscandia_raw.rds")

1.8.1 Distribution of variables

The distribution of the variables were checked with histograms.

fennoscandia_raw <- readRDS("fennoscandia_raw.rds")

h_list <- c()

for (i in names(fennoscandia_raw[,-c(1:2)])){
  if(IQR(fennoscandia_raw[[i]], na.rm = T) > 0.1){
    hi <- ggplot(data = fennoscandia_raw,aes_string(i))+geom_histogram(stat = "bin",na.rm = T,binwidth = function(x) 2*IQR(x, na.rm = T)/(length(x)^(1/3)))+labs(y="",title = i)+theme_light(base_size = 20)
  }else{
    hi <- ggplot(data = fennoscandia_raw,aes_string(i))+geom_histogram(stat = "bin",na.rm = T)+labs(y="",title = i)+theme_light(base_size = 20)
  }
  h_list[[i]] <- hi
}

cowplot::plot_grid(plotlist = h_list,ncol=3)

1.8.2 Transformation of variables

TOC, Runoff, Precipitation and Slope are skewed to the right so they are log-transformed. N-dep, S-dep, Bog and Arable included many zero, so they were not transformed even if they were skewed.

Moreover, all the variables were scaled and centered.

fennoscandia <- readRDS("fennoscandia_raw.rds")

v <- c("toc_mg_l","runoff","temp","prec","slope","alt","bogs","arable","forest","bare","tndep","tsdep")
names(fennoscandia)[which(names(fennoscandia) %in% v)] <- c("TOC", "Runoff","Temp","Precip","Slope","Alt","Bog","Arable","Forest","Bare","TNdep","TSdep")

fennoscandia$logTOC <- log(fennoscandia$TOC) 
fennoscandia$s.logTOC <- scale(fennoscandia$logTOC) 
fennoscandia$Runoff <- fennoscandia$Runoff*365*24*60*60 # from kg/m2/s to kg/m2/year
fennoscandia$logRunoff <- log(fennoscandia$Runoff)  
fennoscandia$s.logRunoff <- scale(fennoscandia$logRunoff)
fennoscandia$NDVI <- fennoscandia$ndvi.value
fennoscandia$s.NDVI <- scale(fennoscandia$NDVI)
fennoscandia$s.Bog <- scale(fennoscandia$Bog)
fennoscandia$s.Arable <- scale(fennoscandia$Arable)
fennoscandia$s.Forest <- scale(fennoscandia$Forest)
fennoscandia$Temp <- fennoscandia$Temp / 10 # divided by 10 because the temperature is stored in C*10
fennoscandia$s.Temp <-  scale(fennoscandia$Temp) 
fennoscandia$logPrecip <- log(fennoscandia$Precip+1e-5)
fennoscandia$s.logPrecip <- scale(fennoscandia$logPrec)
fennoscandia$s.Slope <- scale(fennoscandia$Slope)
fennoscandia$s.TNdep <- scale(fennoscandia$TNdep)
fennoscandia$s.TSdep <- scale(fennoscandia$TSdep)
fennoscandia$rdn <- fennoscandia$drdn+fennoscandia$wrdn
fennoscandia$oxn <- fennoscandia$doxn+fennoscandia$woxn
fennoscandia$s.rdn <- scale(fennoscandia$rdn)
fennoscandia$s.oxn <- scale(fennoscandia$oxn)

# removes old ndvi column
fennoscandia$ndvi <- NULL

saveRDS(fennoscandia, "fennoscandia.rds")

# create fennoscandia SpatialPointsDataFrame
fennoscandia.spdf <- SpatialPointsDataFrame(fennoscandia[,c("longitude","latitude")], fennoscandia)
saveRDS(fennoscandia.spdf, "fennoscandia.spdf.rds")

norge <- filter(fennoscandia,nation == "Norway")
saveRDS(norge, "norge.rds")

The coordinate system of the dataset was converted to UTM33.

catchment.poly <- readRDS("catchment.poly.Rdata")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fennoscandia.cat <- merge(fennoscandia,catchment.poly, by = "ebint") %>% st_as_sf()
fennoscandia.33 <- st_transform(fennoscandia.cat,st_crs(wr.sf.95))
saveRDS(fennoscandia.33, "fennoscandia.33.rds")

2 Mapping of predictor and response variable

The TOC concentration and the predictor variables were plotted by catchment.

fennoscandia.33 <- readRDS("fennoscandia.33.rds")
m <- c("TOC","NDVI","Runoff","Precip","Temp","Bog","Arable","Forest","TNdep","TSdep")
u <- c("TOC concentration,\nmgC/L", "a) NDVI on scale 0 to 1", "b) Mean Runoff, mm/y","c) Mean annual\nprecipitation, mm/y","d) Mean annual\n temperature, °C","e) Proportion of bogs","f) Proportion of\narable land","g) Proportion of forest","h) Total nitrogen\ndeposition, mg/m2", "i) Total sulphur\ndeposition, mg/m2")
couleurs <- c("RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu","RdYlBu")
dir <- c(T,F,T,T,F,F,F,F,F,F)
marg <- 0
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(crs(fennoscandia.33))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(crs(fennoscandia.33))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(crs(fennoscandia.33))

for(i in m[2:10]){
  
 mysf <- fennoscandia.33
 mypal <- brewer.pal(name = couleurs[which(m == i)], n = 9)
 display.brewer.pal(name = couleurs[which(m == i)], n = 9)
 midcol <- median(mysf[[i]])

 if(dir[which(m == i)] == F){
 map <- ggplot()+geom_sf(data = nor, fill = "white", size = 0.2)+
   geom_sf(data = swe, fill = "white", size = 0.2)+
   geom_sf(data = fin, fill = "white", size = 0.2)+
   geom_sf(data = mysf, aes(fill = .data[[i]], col = .data[[i]]))+
   scale_fill_gradientn(colours = rev(mypal),name = paste(u[grep(i,m)],"\n(",i,")",sep=""), aesthetics = c("colour","fill"))+
    theme_minimal(base_size = 8)+
    theme(legend.position = "bottom", plot.margin = unit(c(0,marg,0,marg), units = "cm"))
    
  filename_i <- paste("map",i,"png",sep = ".")
  ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
  
 } else if(dir[which(m == i)] == T){
   map <- ggplot()+geom_sf(data = nor, fill = "white", size = 0.2)+
     geom_sf(data = swe, fill = "white", size = 0.2)+
     geom_sf(data = fin, fill = "white", size = 0.2)+
     geom_sf(data = mysf, aes(fill = .data[[i]], col = .data[[i]]))+
    scale_fill_gradientn(colours = mypal,name = paste(u[grep(i,m)],"\n(",i,")",sep=""), aesthetics = c("colour","fill"))+
                    theme_minimal(base_size = 8)+
    theme(legend.position = "bottom", plot.margin = unit(c(0,marg,0,marg), units = "cm"))
    
  filename_i <- paste("map",i,"png",sep = ".")
  ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
   
 }
}

label_tb <- read.table(text = 
  " Label Long Lat
    Norway 7 65
    Sweden 10 55
    Finland 25 59", header = T) 
label_sf <- label_tb %>% st_as_sf(coords = c("Long", "Lat"), crs = "+proj=longlat +datum=WGS84") %>% st_transform(crs(fennoscandia.33))

for(i in m[1]){
  
 mysf <- fennoscandia.33
 mypal <- brewer.pal(name = couleurs[which(m == i)], n = 9)
 display.brewer.pal(name = couleurs[which(m == i)], n = 9)
 midcol <- quantile(mysf[[i]], 0.99)

 map <- ggplot()+geom_sf(data = nor, fill = "white", size = 0.2)+
  geom_sf_label(data = label_sf, aes(label = Label), label.size = 0.2, label.padding = unit(0.1,"lines"), size = 3)+
  geom_sf(data = swe, fill = "white", size = 0.2)+
  geom_sf(data = fin, fill = "white", size = 0.2)+
  geom_sf(data = filter(mysf, TOC < 20), aes(fill = .data[[i]], col = .data[[i]]))+
  geom_sf(data = filter(mysf, TOC > 20), fill = "sienna4", col = "sienna4")+
  scale_fill_gradientn(colours = rev(mypal),name = paste(u[grep(i,m)],"\n(",i,")",sep=""), aesthetics = c("colour","fill"))+
   labs(x="",y="")+
   theme_minimal(base_size = 9)+
    theme(legend.position = "bottom", plot.margin = unit(c(0,marg,0,marg), units = "cm"))
    
  filename_i <- paste("map",i,"png",sep = ".")
  ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
  
}

white.fen <- ggplot()+theme_void()+theme(plot.margin = unit(c(0,marg,0,marg), units = "cm"))
ggsave(plot = white.fen, filename = "white.fen.png", dpi = "retina", width = 10, height = 12, units = "cm")
knitr::include_graphics(c("map.TOC.png","map.NDVI.png","map.Runoff.png","map.Precip.png","map.Temp.png","map.Bog.png","map.Arable.png","map.Forest.png","map.TNdep.png","map.TSdep.png"))

3 Spatial autocorrelation & VIF

The spatial autocorrelation of the response and predictor variables were computed using Moran’s I (see Supplementary 4).

The neighbor matrix is computed first.

fennoscandia.spdf <- readRDS("fennoscandia.spdf.rds")
k.neigh.w <- knearneigh(fennoscandia.spdf, k = 100) %>% knn2nb() %>% nb2listw()
saveRDS(k.neigh.w,"k.neigh.w.rds")

The spatial autocorrelation is computed for the response and predictor variables.

k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")

moran.list <- c() 

m <- c("logTOC","NDVI","logRunoff","logPrecip","Temp","Bog","Arable","Forest","TNdep","TSdep")

for (i in m){
  moran.list[[i]] <- moran.test(fennoscandia[,i],k.neigh.w, na.action = na.omit, alternative = "two.sided")
}

moran.df <- data.frame(m) %>% setNames("parameter")
moran.df$moran.I <- NA
moran.df$p.value <- NA

for (i in 1:length(moran.list)){
  moran.df$moran.I[i] <- moran.list[[i]]$estimate[1]
  moran.df$p.value[i] <- moran.list[[i]]$p.value
}

saveRDS(moran.df,"moran.df.rds")
moran.df <- readRDS("moran.df.rds")
moran.df %>% arrange(desc(moran.I)) %>% kable %>% kable_styling(bootstrap_options = "bordered", position = "center")
parameter moran.I p.value
Temp 0.8975763 0
logPrecip 0.8893982 0
TNdep 0.8781398 0
TSdep 0.8580784 0
logRunoff 0.7907892 0
logTOC 0.5796447 0
Forest 0.4847911 0
NDVI 0.4622196 0
Bog 0.1741102 0
Arable 0.1275870 0

4 Comparison LM and SELM, 5 selected predictors (NDVI, Runoff, Bog, Arable, TNdep)

4.1 Linear model on NSF, 5 selected predictors

A linear model was fitted with the 5 selected predictors, both for scaled and unscaled variables.

k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")

fm.s <- s.logTOC~s.NDVI+s.logRunoff+s.Bog+s.Arable+s.TNdep
s.lm.fennoscandia.5 <- lm(fm.s, data = fennoscandia, na.action = na.exclude)
saveRDS(s.lm.fennoscandia.5,"s.lm.fennoscandia.5.rds")

fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
lm.fennoscandia.5 <- lm(fm, data = fennoscandia, na.action = na.exclude)
saveRDS(lm.fennoscandia.5,"lm.fennoscandia.5.rds")

s.lm.r2.5 <- summary(s.lm.fennoscandia.5)$adj.r.squared
s.moran.I.res.lm.5 <- lm.morantest(s.lm.fennoscandia.5, k.neigh.w, alternative = "two.sided")

lm.r2.5 <- summary(lm.fennoscandia.5)$adj.r.squared
moran.I.res.lm.5 <- lm.morantest(lm.fennoscandia.5, k.neigh.w, alternative = "two.sided")

moran.limit.fennoscandia.5 <- 1/(dim(fennoscandia)[1]-1)

s.lm.coef.5 <- summary(s.lm.fennoscandia.5)$coefficients %>% as.data.frame()
s.lm.coef.5 <- s.lm.coef.5[-1,] %>% rownames_to_column(var = "Parameter")
saveRDS(s.lm.coef.5, "s.lm.coef.5.rds")

lm.coef.5 <- summary(lm.fennoscandia.5)$coefficients %>% as.data.frame()
lm.coef.5 <- lm.coef.5[-1,] %>% rownames_to_column(var = "Parameter")

lm.effectsize.5 <- compute.effect.size(lm.coef.5)
saveRDS(lm.effectsize.5, "lm.effectsize.5.rds")

We checked the distribution of the coefficients, a posteriori, using SIM.

sm.fennoscandia.5 <- arm::sim(lm.fennoscandia.5,n.sims = 100)

sm.coef.5 <- sm.fennoscandia.5@coef[,-1]
sm.effectsize.5 <- data.frame(parameter = colnames(sm.coef.5), effect.size = NA, std.error = NA)

n <- dim(sm.coef.5)[2]
m <- colnames(sm.coef.5)

for(i in 1:n){
  if(grepl("log",m[i])==T){
    sm.effectsize.5[i,2] <- mean(1.1^(sm.coef.5[,i]))
    sm.effectsize.5[i,3] <- sd(1.1^(sm.coef.5[,i])*log(1.1)*sm.coef.5[,i])
  }else if(m[i] %in% c("NDVI","Bog","Arable","Forest")){
    sm.effectsize.5[i,2] <- mean(exp(sm.coef.5[,i]*0.1))
    sm.effectsize.5[i,3] <- sd(0.1*exp(sm.coef.5[,i]*0.1)*sm.coef.5[,i])
  }else if(m[i] %in% c("TNdep","TSdep")){
    sm.effectsize.5[i,2] <- mean(exp(sm.coef.5[,i]*250))
    sm.effectsize.5[i,3] <- sd(250*exp(sm.coef.5[,i]*250)*sm.coef.5[,i])
  }
}

The results are summarized below.

summary.table.lm.5 <- cbind(dplyr::select(lm.coef.5,Parameter),s.lm.coef.5[,2:3],lm.coef.5[,2:3],lm.effectsize.5[,2:4],sm.effectsize.5[,2:3]) %>%
                rbind(c("AIC",AIC(s.lm.fennoscandia.5),"",AIC(lm.fennoscandia.5),"","","","","","")) %>%
                rbind(c("R2",s.lm.r2.5,"",lm.r2.5,"","","","","","")) %>%
                rbind(c("Moran'I res",s.moran.I.res.lm.5$estimate[[1]],"",moran.I.res.lm.5$estimate[[1]],"","","","","",""))
summary.table.lm.5[,-1] <- apply(summary.table.lm.5[,-1],2,as.numeric)

saveRDS(summary.table.lm.5, "summary.table.lm.5.rds")
summary.table.lm.5 <- readRDS("summary.table.lm.5.rds")

n <- 5

knitr::kable(summary.table.lm.5,digits = 3, caption = "LM with 5 selected predictors") %>%
  add_header_above(c("Model"=1,"Scaled LM" = 2, "LM" = 2, "Effect size" = 3,"SIM" = 2),italic = T) %>%
  column_spec(column = 8, bold = T) %>% 
  kable_styling(bootstrap_options = "bordered") %>%
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+3)
LM with 5 selected predictors
Model
Scaled LM
LM
Effect size
SIM
Parameter Estimate Std. Error Estimate Std. Error Estimate Std. Error Percent effect.size std.error
NDVI 0.442 0.010 3.355 0.079 1.034 0.003 3.412 1.400 0.015
logRunoff -0.429 0.011 -0.762 0.019 0.992 0.000 -0.755 0.930 0.001
Bog 0.124 0.009 0.799 0.059 1.008 0.000 0.803 1.083 0.007
Arable 0.002 0.009 0.027 0.123 1.000 0.000 0.027 1.005 0.013
TNdep 0.082 0.010 0.000 0.000 1.004 0.001 0.445 1.046 0.006
Model Evaluation
AIC 8809.136 7241.057
R2 0.624 0.624
Moran’I res 0.154 0.154

rtable_nums(“summary-lm-5”,“Results for the linear model with 5 selected predictors on Norway, Sweden and Finland”)`

The residuals are plotted below.

s.lm.fennoscandia.5 <- readRDS("s.lm.fennoscandia.5.rds")
lm.fennoscandia.5 <- readRDS("lm.fennoscandia.5.rds")

f_plotspatial(data = fennoscandia,var = s.lm.fennoscandia.5$residuals, plottitle = "a) Residuals for scaled LM")

f_plotspatial(data = fennoscandia,var = lm.fennoscandia.5$residuals, plottitle = "b) Residuals for unscaled LM")

rfig_nums(“residuals-lm-5”,“Residuals for linear model, a) with scaled variables, b) with unscaled variables”)`

4.2 Spatial error linear model on NSF, 5 selected predictors (NDVI, Runoff, Bog, Arable, TNdep)

A Spatial Error Linear Model was fitted on for the catchments in Norway, Sweden and Finland, both on scaled and unscaled variables.

k.neigh.w <- readRDS("k.neigh.w.rds")
norge.kmat <- readRDS("norge.kmat.Rdata")

#lagrange <- lm.LMtests(lm.fennoscandia,k.neigh.w,test = c("LMerr","LMlag","RLMerr","RLMlag"))

fm.s <- s.logTOC~s.NDVI+s.logRunoff+s.Bog+s.Arable+s.TNdep
s.sem.fennoscandia.5 <- errorsarlm(fm.s,fennoscandia,k.neigh.w)
saveRDS(s.sem.fennoscandia.5,"s.sem.fennoscandia.5.rds")

fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
sem.fennoscandia.5 <- errorsarlm(fm,fennoscandia,k.neigh.w)
saveRDS(sem.fennoscandia.5,"sem.fennoscandia.5.rds")

The results of the model are displayed below.

s.sem.fennoscandia.5 <- readRDS("s.sem.fennoscandia.5.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.rds")
k.neigh.w <- readRDS("k.neigh.w.rds")

s.sem.coef.5 <- s.sem.fennoscandia.5$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
s.sem.coef.5 <- s.sem.coef.5 %>% rownames_to_column(var = "Parameter")
s.sem.coef.5$std.error <- s.sem.fennoscandia.5$rest.se[-1]
saveRDS(s.sem.coef.5, "s.sem.coef.5.rds")

sem.coef.5 <- sem.fennoscandia.5$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
sem.coef.5 <- sem.coef.5 %>% rownames_to_column(var = "Parameter")
sem.coef.5$std.error <- sem.fennoscandia.5$rest.se[-1]


sem.effectsize.5 <- compute.effect.size(sem.coef.5)
saveRDS(sem.effectsize.5, "sem.effectsize.5.rds")

n <- length(sem.coef.5$Parameter)

s.sem.AIC.5 <- 2*(n-1)-2*s.sem.fennoscandia.5$LL
sem.AIC.5 <- 2*(n-1)-2*sem.fennoscandia.5$LL

s.moran.I.res.sem.5 <- moran.test(s.sem.fennoscandia.5$residuals, k.neigh.w, alternative = "two.sided")

moran.I.res.sem.5 <- moran.test(sem.fennoscandia.5$residuals, k.neigh.w, alternative = "two.sided")

summary.table.sem.5 <- cbind(dplyr::select(sem.coef.5,Parameter),s.sem.coef.5[,2:3],sem.coef.5[,2:3],sem.effectsize.5[,2:4]) %>%
                rbind(c("AIC",s.sem.AIC.5,"",sem.AIC.5,"","","","")) %>%
                rbind(c("Moran'I res", s.moran.I.res.sem.5$estimate[[1]],"", moran.I.res.sem.5$estimate[[1]], "","","",""))
summary.table.sem.5[,-1] <- apply(summary.table.sem.5[,-1],2,as.numeric)

saveRDS(summary.table.sem.5, "summary.table.sem.5.rds")
summary.table.se.5 <- readRDS("summary.table.sem.5.rds")

n <- 5

knitr::kable(summary.table.sem.5,digits = 3, caption = "SELM with 5 selected predictors") %>%
  add_header_above(c("Model"=1,"Scaled SELM" = 2, "SELM" = 2, "Effect size" = 3),italic = T) %>%
  kable_styling(bootstrap_options = "bordered") %>%
  column_spec(column = 8, bold = T) %>% 
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+2)
SELM with 5 selected predictors
Model
Scaled SELM
SELM
Effect size
Parameter Estimate std.error Estimate std.error Estimate std.error Percent
NDVI 0.327 0.012 2.480 0.088 1.025 0.002 2.511
logRunoff -0.181 0.022 -0.322 0.038 0.997 0.000 -0.320
Bog 0.144 0.009 0.932 0.056 1.009 0.001 0.936
Arable -0.003 0.009 -0.041 0.114 1.000 0.000 -0.041
TNdep 0.125 0.031 0.000 0.000 1.007 0.002 0.678
Model Evaluation
AIC 7742.975 6174.896
Moran’I res 0.010 0.010

The residuals for scaled and unscaled models are presented in Figure 1

s.sem.fennoscandia.5 <- readRDS("s.sem.fennoscandia.5.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.rds")

f_plotspatial(data = fennoscandia,var = s.sem.fennoscandia.5$residuals, plottitle = "a) Residuals for scaled SELM")

f_plotspatial(data = fennoscandia,var = sem.fennoscandia.5$residuals, plottitle = "b)Residuals for unscaled SELM")

Figure 1: Map of residuals for a) scaled and b) unscaled SELM

4.3 Comparative plots 5

s.lm.coef.5 <- readRDS("s.lm.coef.5.rds")
s.sem.coef.5 <- readRDS("s.sem.coef.5.rds")

n <- c(names(s.sem.coef.5),"Model")
m <- cbind.data.frame(s.lm.coef.5[,1:3],Model = "LM") %>% setNames(n)
o <- cbind.data.frame(s.sem.coef.5,Model = "SELM")
  
scaled.df.5 <- rbind(m,o)
scaled.df.5$Parameter <- factor(scaled.df.5$Parameter, levels = c("s.TNdep","s.Arable","s.Bog","s.logRunoff","s.NDVI"))

scaled.plot.5 <- ggplot(scaled.df.5, aes(group = Model))+geom_col(aes(x=Estimate,y=Parameter, fill = Model), position = "dodge", width = 0.5)+
  scale_fill_manual(values = c("#d73027","#4575b4"))+
  geom_errorbar(aes(y=Parameter,xmin = Estimate-std.error,xmax=Estimate+std.error), position = "dodge", width = 0.5)+
  labs(x = "Scaled Estimate", y = "Model with 6 predictors")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
sem.effectsize.5 <- readRDS("sem.effectsize.5.rds")
lm.effectsize.5 <- readRDS("lm.effectsize.5.rds")

n <- c(names(sem.effectsize.5),"Model")
m <- cbind.data.frame(lm.effectsize.5,Model = "LM") %>% setNames(n)
o <- cbind.data.frame(sem.effectsize.5,Model = "SELM")
  
h.df.5 <- rbind(m,o)
h.df.5$Parameter <- factor(h.df.5$Parameter, levels = c("TNdep","Arable","Bog","logRunoff","NDVI"))

effectsize.plot.5 <- ggplot(h.df.5, aes(group = Model))+
  geom_col(aes(x=Percent,y=Parameter, fill = Model), position = "dodge", width = 0.5)+
  scale_fill_manual(values = c("#d73027","#4575b4"))+
  geom_errorbar(aes(y=Parameter,xmin = Percent-(100*std.error),xmax=Percent+(100*std.error)), position = "dodge", width = 0.5)+
  labs(x = "Effect size in %", y = "")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
compare.5 <- cowplot::plot_grid(plotlist = list(scaled.plot.5,effectsize.plot.5),ncol=2)
save_plot(plot = compare.5, filename = "compare.lm.selm.5.png")
print(compare.5)

Figure 2: Comparison of scaled estimates for LM and SELM (left) and of the effect sizes for each predictor (right) with NDVI, logRunoff, Bog, Arable and TNdep as predictors

5 Other models

5.1 Comparison of linear model and spatial error model, 9 predictors

The linear model and the spatial error linear model were fitted with all 9 predictors (including predictors correlated to each other or having a high Variance Inflation Factor).

5.1.1 Linear model on NSF, 9 predictors

k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")

fm.s <- s.logTOC~s.NDVI+s.logRunoff+s.Bog+s.Arable+s.Forest+s.TNdep+s.TSdep+s.Temp+s.logPrecip
s.lm.fennoscandia.9 <- lm(fm.s, data = fennoscandia, na.action = na.exclude)
saveRDS(s.lm.fennoscandia.9,"s.lm.fennoscandia.9.rds")

fm <- logTOC~NDVI+logRunoff+Bog+Arable+Forest+TNdep+TSdep+Temp+logPrecip
lm.fennoscandia.9 <- lm(fm, data = fennoscandia, na.action = na.exclude)
saveRDS(lm.fennoscandia.9,"lm.fennoscandia.9.rds")

s.lm.r2 <- summary(s.lm.fennoscandia.9)$adj.r.squared
s.moran.I.res.lm <- lm.morantest(s.lm.fennoscandia.9, k.neigh.w, alternative = "two.sided")

lm.r2 <- summary(lm.fennoscandia.9)$adj.r.squared
moran.I.res.lm <- lm.morantest(lm.fennoscandia.9, k.neigh.w, alternative = "two.sided")

moran.limit.fennoscandia <- 1/(dim(fennoscandia)[1]-1)

s.lm.coef.9 <- summary(s.lm.fennoscandia.9)$coefficients %>% as.data.frame()
s.lm.coef.9 <- s.lm.coef.9[-1,] %>% rownames_to_column(var = "Parameter")
saveRDS(s.lm.coef.9, "s.lm.coef.9.rds")

lm.coef.9 <- summary(lm.fennoscandia.9)$coefficients %>% as.data.frame()
lm.coef.9 <- lm.coef.9[-1,] %>% rownames_to_column(var = "Parameter")

lm.effectsize.9 <- compute.effect.size(lm.coef.9)
saveRDS(lm.effectsize.9, "lm.effectsize.9.rds")
summary.table.lm.9 <- cbind(dplyr::select(lm.coef.9,Parameter),s.lm.coef.9[,2:3],lm.coef.9[,2:3],lm.effectsize.9[,2:4]) %>%
                rbind(c("AIC",AIC(s.lm.fennoscandia.9),"",AIC(lm.fennoscandia.9),"","","")) %>%
                rbind(c("R2",s.lm.r2,"",lm.r2,"","","")) %>%
                rbind(c("Moran'I res",s.moran.I.res.lm$estimate[[1]],"",moran.I.res.lm$estimate[[1]],"","",""))
summary.table.lm.9[,-1] <- apply(summary.table.lm.9[,-1],2,as.numeric)

saveRDS(summary.table.lm.9, "summary.table.lm.9.rds")

Table 1: Summary table for linear model with 9 predictors

summary.table.lm.9 <- readRDS("summary.table.lm.9.rds")

n <- 9

knitr::kable(summary.table.lm.9,digits = 3, caption = "LM with 9 predictors") %>%
  add_header_above(c("Model"=1,"Scaled LM" = 2, "LM" = 2, "Effect size" = 3),italic = T) %>%
  column_spec(column = 8, bold = T) %>% 
  kable_styling(bootstrap_options = "bordered") %>%
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+3)
LM with 9 predictors
Model
Scaled LM
LM
Effect size
Parameter Estimate Std. Error Estimate Std. Error Estimate Std. Error Percent
NDVI 0.278 0.011 2.109 0.081 1.021 0.002 2.132
logRunoff -0.006 0.018 -0.011 0.033 1.000 0.000 -0.011
Bog 0.175 0.009 1.133 0.057 1.011 0.001 1.140
Arable 0.021 0.009 0.270 0.117 1.003 0.000 0.271
Forest 0.161 0.013 0.448 0.036 1.004 0.000 0.449
TNdep -0.186 0.033 0.000 0.000 0.990 0.002 -1.000
TSdep 0.122 0.036 0.001 0.000 1.007 0.002 0.677
Temp 0.449 0.019 0.141 0.006 1.001 0.000 0.141
logPrecip -0.395 0.018 -1.097 0.050 0.989 0.000 -1.085
Model Evaluation
AIC 7666.819 6098.740
R2 0.705 0.705
Moran’I res 0.101 0.101 
s.lm.fennoscandia.9 <- readRDS("s.lm.fennoscandia.9.rds")
lm.fennoscandia.9 <- readRDS("lm.fennoscandia.9.rds")

f_plotspatial(data = fennoscandia,var = s.lm.fennoscandia.9$residuals, plottitle = "Residuals for scaled LM, with all predictors")

f_plotspatial(data = fennoscandia,var = lm.fennoscandia.9$residuals, plottitle = "Residuals for unscaled LM, with all predictors")

Figure 3: Map of residuals for a) scaled and b) unscaled LM with 9 predictors

5.1.2 Spatial error linear model on NSF, 9 predictors

k.neigh.w <- readRDS("k.neigh.w.rds")

fm.s <- s.logTOC~s.NDVI+s.logRunoff+s.Bog+s.Arable+s.Forest+s.TNdep+s.TSdep+s.Temp+s.logPrecip
s.sem.fen.9 <- errorsarlm(fm.s,fennoscandia,k.neigh.w)
saveRDS(s.sem.fen.9,"s.sem.fen.9.Rdata")

fm <- logTOC~NDVI+logRunoff+Bog+Arable+Forest+TNdep+TSdep+Temp+logPrecip
sem.fen.9 <- errorsarlm(fm,fennoscandia,k.neigh.w)
saveRDS(sem.fen.9,"sem.fen.9.Rdata")
s.sem.fen.9 <- readRDS("s.sem.fen.9.Rdata")
sem.fen.9 <- readRDS("sem.fen.9.Rdata")
k.neigh.w <- readRDS("k.neigh.w.rds")


s.sem.coef <- s.sem.fen.9$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
s.sem.coef <- s.sem.coef %>% rownames_to_column(var = "Parameter")
s.sem.coef$std.error <- s.sem.fen.9$rest.se[-1]
saveRDS(s.sem.coef, "s.sem.coef.9.rds")

sem.coef <- sem.fen.9$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
sem.coef <- sem.coef %>% rownames_to_column(var = "Parameter")
sem.coef$std.error <- sem.fen.9$rest.se[-1]


sem.effectsize <- compute.effect.size(sem.coef)
saveRDS(sem.effectsize, "sem.effectsize.9.rds")

n <- length(sem.coef$Parameter)

s.sem.AIC <- 2*(n-1)-2*s.sem.fen.9$LL
sem.AIC <- 2*(n-1)-2*sem.fen.9$LL

s.moran.I.res.sem <- moran.test(s.sem.fen.9$residuals, k.neigh.w, alternative = "two.sided")

moran.I.res.sem <- moran.test(sem.fen.9$residuals, k.neigh.w, alternative = "two.sided")

summary.table.sem <- cbind(dplyr::select(sem.coef,Parameter),s.sem.coef[,2:3],sem.coef[,2:3],sem.effectsize[,2:4]) %>%
                rbind(c("AIC",s.sem.AIC,"",sem.AIC,"","","","")) %>%
                rbind(c("Moran'I res", s.moran.I.res.sem$estimate[[1]],"", moran.I.res.sem$estimate[[1]], "","","",""))
summary.table.sem[,-1] <- apply(summary.table.sem[,-1],2,as.numeric)

saveRDS(summary.table.sem, "summary.table.sem.9.rds")

Table 2: Summary table for spatial error linear model with 9 predictors

summary.table.sem.9 <- readRDS("summary.table.sem.9.rds")

n <- 9

knitr::kable(summary.table.sem.9,digits = 3, caption = "SELM with 9 predictors") %>%
  add_header_above(c("Model"=1,"Scaled SELM" = 2, "SELM" = 2, "Effect size" = 3),italic = T) %>%
  kable_styling(bootstrap_options = "bordered") %>%
  column_spec(column = 8, bold = T) %>% 
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+2)
SELM with 9 predictors
Model
Scaled SELM
SELM
Effect size
Parameter Estimate std.error Estimate std.error Estimate std.error Percent
NDVI 0.189 0.012 1.434 0.091 1.014 0.001 1.445
logRunoff -0.034 0.022 -0.060 0.039 0.999 0.000 -0.059
Bog 0.163 0.009 1.052 0.056 1.011 0.001 1.058
Arable 0.000 0.009 -0.002 0.117 1.000 0.000 -0.002
Forest 0.139 0.013 0.387 0.035 1.004 0.000 0.388
TNdep -0.126 0.073 0.000 0.000 0.993 0.004 -0.679
TSdep -0.055 0.071 0.000 0.000 0.997 0.004 -0.305
Temp 0.633 0.036 0.198 0.011 1.002 0.000 0.199
logPrecip -0.425 0.032 -1.182 0.090 0.988 0.001 -1.169
Model Evaluation
AIC 7109.149 5541.070
Moran’I res 0.009 0.009 
s.sem.fen.9 <- readRDS("s.sem.fen.9.Rdata")
sem.fen.9 <- readRDS("sem.fen.9.Rdata")

f_plotspatial(data = fennoscandia,var = s.sem.fen.9$residuals, plottitle = "a) Residuals for scaled SEM, with all predictors")

f_plotspatial(data = fennoscandia,var = sem.fen.9$residuals, plottitle = "b) Residuals for unscaled SELM, with all predictors")

Figure 4: Map of residuals for a) scaled and b) unscaled SELM with 9 predictors

5.1.3 Comparative plots 9

s.lm.coef.9 <- readRDS("s.lm.coef.9.rds")
s.sem.coef.9 <- readRDS("s.sem.coef.9.rds")

n <- c(names(s.sem.coef.9),"Model")
m <- cbind.data.frame(s.lm.coef.9[,1:3],Model = "LM") %>% setNames(n)
o <- cbind.data.frame(s.sem.coef.9,Model = "SELM")
  
scaled.df <- rbind(m,o)
scaled.df$Parameter <- factor(scaled.df$Parameter, levels = c("s.TNdep","s.TSdep","s.Forest","s.Arable","s.Bog","s.logRunoff","s.Temp","s.logPrecip","s.NDVI"))

scaled.plot.9 <- ggplot(scaled.df, aes(group = Model))+geom_col(aes(x=Estimate,y=Parameter, fill = Model), position = "dodge", width = 0.5)+ scale_fill_manual(values = c("#d73027","#4575b4"))+
  geom_errorbar(aes(y=Parameter,xmin = Estimate-std.error,xmax=Estimate+std.error), position = "dodge", width = 0.5)+
  labs(x = "Scaled Estimate", y = "Model with 9 predictors")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
lm.effectsize.9 <- readRDS("lm.effectsize.9.rds")
sem.effectsize.9 <- readRDS("sem.effectsize.9.rds")

n <- c(names(sem.effectsize.9),"Model")
m <- cbind.data.frame(lm.effectsize.9,Model = "LM") %>% setNames(n)
o <- cbind.data.frame(sem.effectsize.9,Model = "SELM")
  
h.df <- rbind(m,o)
h.df$Parameter <- factor(h.df$Parameter, levels = c("TNdep","TSdep","Forest","Arable","Bog","logRunoff","Temp","logPrecip","NDVI"))

effectsize.plot.9 <- ggplot(h.df, aes(group = Model))+
  geom_col(aes(x=Percent,y=Parameter, fill = Model), position = "dodge", width = 0.5)+ scale_fill_manual(values = c("#d73027","#4575b4"))+
  geom_errorbar(aes(y=Parameter,xmin = Percent-(100*std.error),xmax=Percent+(100*std.error)), position = "dodge", width = 0.5)+
  labs(x = "Effect size in %", y = "")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
compare.9 <- cowplot::plot_grid(plotlist = list(scaled.plot.9,effectsize.plot.9),ncol=2)
save_plot(plot = compare.9, filename = "compare.lm.selm.9.png")
print(compare.9)

Figure 5: Comparison of scaled estimates for LM and SELM (left) and of the effect sizes for each predictor (right) with 9 predictors

5.2 Spatial error linear models with NDVI, Bog and logRunoff

5.3 Spatial error linear model on Fennoscandian datasets with NDVI, Bog and logRunoff

k.neigh.w <- readRDS("k.neigh.w.Rdata")
fennoscandia <- readRDS("fennoscandia.rds")

s.fm <- s.logTOC~s.NDVI+s.logRunoff+s.Bog
s.sem.fen.3 <- errorsarlm(s.fm,fennoscandia,k.neigh.w)
saveRDS(s.sem.fen.3,"s.sem.fen.3.rds")

fm <- logTOC~NDVI+logRunoff+Bog
sem.fen.3 <- errorsarlm(fm,fennoscandia,k.neigh.w)
saveRDS(sem.fen.3,"sem.fen.3.rds")

s.sem.fen.3 <- readRDS("s.sem.fen.3.rds")
sem.fen.3 <- readRDS("sem.fen.3.rds")
k.neigh.w <- readRDS("k.neigh.w.Rdata")

s.sem.coef.3 <- s.sem.fen.3$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
s.sem.coef.3 <- s.sem.coef.3 %>% rownames_to_column(var = "Parameter")
s.sem.coef.3$std.error <- s.sem.fen.3$rest.se[-1]
saveRDS(s.sem.coef.3, "s.sem.coef.3.rds")

sem.coef.3 <- sem.fen.3$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
sem.coef.3 <- sem.coef.3 %>% rownames_to_column(var = "Parameter")
sem.coef.3$std.error <- sem.fen.3$rest.se[-1]

sem.effectsize.fen.3 <- compute.effect.size(sem.coef.3)
saveRDS(sem.effectsize.fen.3, "sem.effectsize.fen.3.rds")


n <- length(sem.coef.3$Parameter)

s.sem.AIC.3 <- 2*(n-1)-2*s.sem.fen.3$LL
s.moran.I.res.sem.3 <- moran.test(s.sem.fen.3$residuals, k.neigh.w, alternative = "two.sided")

sem.AIC.3 <- 2*(n-1)-2*sem.fen.3$LL
moran.I.res.sem.3 <- moran.test(sem.fen.3$residuals, k.neigh.w, alternative = "two.sided")

summary.table.simple.sem.3 <- cbind(dplyr::select(sem.coef.3,Parameter),s.sem.coef.3[,2:3],sem.coef.3[,2:3],sem.effectsize.fen.3[,2:4]) %>%
                rbind(c("AIC",s.sem.AIC.3,"",sem.AIC.3,"","","","")) %>%
                rbind(c("Moran'I res", s.moran.I.res.sem.3$estimate[[1]],"",moran.I.res.sem.3$estimate[[1]], "","","",""))
summary.table.simple.sem.3[,-1] <- apply(summary.table.simple.sem.3[,-1],2,as.numeric)

saveRDS(summary.table.simple.sem.3, "summary.table.simple.sem.3.rds")

Table 3: Summary table for SELM model with 3 predictors

summary.table.simple.sem.3 <- readRDS("summary.table.simple.sem.3.rds")
n <- 3

knitr::kable(summary.table.simple.sem.3,digits = 3, caption = "SELM with 3 predictors (NDVI, Bog, Runoff") %>%
  add_header_above(c("Model"=1,"Scaled SELM" = 2,"SELM" = 2, "Effect size" = 3),italic = T) %>%
  kable_styling(bootstrap_options = "bordered") %>%
  column_spec(column = 8, bold = T) %>% 
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+2)
SELM with 3 predictors (NDVI, Bog, Runoff
Model
Scaled SELM
SELM
Effect size
Parameter Estimate std.error Estimate std.error Estimate std.error Percent
NDVI 0.331 0.012 2.510 0.087 1.025 0.002 2.542
logRunoff -0.220 0.019 -0.391 0.034 0.996 0.000 -0.388
Bog 0.143 0.009 0.927 0.057 1.009 0.001 0.931
Model Evaluation
AIC 7755.641 6187.561
Moran’I res 0.011 0.011 
s.sem.fennoscandia.3 <- readRDS("s.sem.fen.3.rds")
sem.fennoscandia.3 <- readRDS("sem.fen.3.rds")

f_plotspatial(data = fennoscandia,var = s.sem.fennoscandia.3$residuals, plottitle = "Residuals for scaled SEM, with NDVI, Runoff and Bogs")

f_plotspatial(data = fennoscandia,var = sem.fennoscandia.3$residuals, plottitle = "Residuals for unscaled SELM, with NDVI, Runoff and Bogs")

Figure 6: Map of residuals for a) scaled and b) unscaled SELM with 3 predictors

sem.coef.3 <- readRDS("s.sem.coef.3.rds")

scaled.plot.3 <- ggplot(s.sem.coef.3)+geom_col(aes(x=Estimate,y=Parameter), width = 0.5, fill = "#d73027")+
  geom_errorbar(aes(y=Parameter,xmin = Estimate-std.error,xmax=Estimate+std.error), width = 0.5)+
  labs(x = "Scaled estimates",y = "Predictor")+
  theme_bw(base_size = 8)+
  theme(legend.position = "bottom")
sem.effectsize.fen.3 <- readRDS("sem.effectsize.fen.3.rds")

effectsize.plot.3 <- ggplot(sem.effectsize.fen.3)+
  geom_col(aes(x=Percent,y=Parameter), width = 0.5, fill = "#4575b4")+
  geom_errorbar(aes(y=Parameter,xmin = Percent-(100*std.error),xmax=Percent+(100*std.error)), position = "dodge", width = 0.5)+
  labs(x = "Effect size in %", y = "")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
compare.3 <- cowplot::plot_grid(plotlist = list(scaled.plot.3,effectsize.plot.3),ncol=2)
save_plot(plot = compare.3, filename = "compare.lm.selm.3.png" )
print(compare.3)

Figure 7: Comparison of scaled estimates for LM and SELM (left) and of the effect sizes for each predictor (right) with 9 predictors

5.4 Spatial error linear model on Fennoscandian datasets with NDVI, Bog and logPrecip

k.neigh.w <- readRDS("k.neigh.w.Rdata")
fennoscandia <- readRDS("fennoscandia.rds")

s.fm <- s.logTOC~s.NDVI+s.logPrecip+s.Bog
s.sem.fen.precip <- errorsarlm(s.fm,fennoscandia,k.neigh.w)
saveRDS(s.sem.fen.precip,"s.sem.fen.precip.rds")

fm <- logTOC~NDVI+logPrecip+Bog
sem.fen.precip <- errorsarlm(fm,fennoscandia,k.neigh.w)
saveRDS(sem.fen.precip,"sem.fen.precip.rds")

s.sem.fen.precip <- readRDS("s.sem.fen.precip.rds")
sem.fen.precip <- readRDS("sem.fen.precip.rds")
k.neigh.w <- readRDS("k.neigh.w.Rdata")

s.sem.coef.precip <- s.sem.fen.precip$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
s.sem.coef.precip <- s.sem.coef.precip %>% rownames_to_column(var = "Parameter")
s.sem.coef.precip$std.error <- s.sem.fen.precip$rest.se[-1]
saveRDS(s.sem.coef.precip, "s.sem.coef.precip.rds")

sem.coef.precip <- sem.fen.precip$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
sem.coef.precip <- sem.coef.precip %>% rownames_to_column(var = "Parameter")
sem.coef.precip$std.error <- sem.fen.precip$rest.se[-1]

sem.effectsize.precip <- compute.effect.size(sem.coef.precip)
saveRDS(sem.effectsize.precip, "sem.effectsize.precip.rds")


n <- length(sem.coef.precip$Parameter)

s.sem.AIC.precip <- 2*(n-1)-2*s.sem.fen.precip$LL
s.moran.I.res.sem.precip <- moran.test(s.sem.fen.precip$residuals, k.neigh.w, alternative = "two.sided")

sem.AIC.precip <- 2*(n-1)-2*sem.fen.precip$LL
moran.I.res.sem.precip <- moran.test(sem.fen.precip$residuals, k.neigh.w, alternative = "two.sided")

summary.table.simple.sem.precip <- cbind(dplyr::select(sem.coef.precip,Parameter),s.sem.coef.precip[,2:3],sem.coef.precip[,2:3],sem.effectsize.precip[,2:4]) %>%
                rbind(c("AIC",s.sem.AIC.precip,"",sem.AIC.precip,"","","","")) %>%
                rbind(c("Moran'I res", s.moran.I.res.sem.precip$estimate[[1]],"",moran.I.res.sem.precip$estimate[[1]], "","","",""))
summary.table.simple.sem.precip[,-1] <- apply(summary.table.simple.sem.precip[,-1],2,as.numeric)

saveRDS(summary.table.simple.sem.precip, "summary.table.simple.precip.rds")
summary.table.simple.precip <- readRDS("summary.table.simple.precip.rds")

n <- 3

knitr::kable(summary.table.simple.precip,digits = 3, caption = "SELM with 3 predictors (NDVI, Bog, logPrecip") %>%
  add_header_above(c("Model"=1,"Scaled SELM" = 2,"SELM" = 2, "Effect size" = 3),italic = T) %>%
  kable_styling(bootstrap_options = "bordered") %>%
  column_spec(column = 8, bold = T) %>% 
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+2)
SELM with 3 predictors (NDVI, Bog, logPrecip
Model
Scaled SELM
SELM
Effect size
Parameter Estimate std.error Estimate std.error Estimate std.error Percent
NDVI 0.323 0.012 2.450 0.088 1.025 0.002 2.480
logPrecip -0.417 0.033 -1.159 0.091 0.989 0.001 -1.147
Bog 0.139 0.009 0.900 0.056 1.009 0.001 0.904
Model Evaluation
AIC 7724.981 6156.902
Moran’I res 0.014 0.014 
s.sem.fen.precip <- readRDS("s.sem.fen.precip.rds")
sem.fen.precip <- readRDS("sem.fen.precip.rds")

f_plotspatial(data = fennoscandia,var = s.sem.fen.precip$residuals, plottitle = "Residuals for scaled SEM, with NDVI, logPrecip, Bogs")

f_plotspatial(data = fennoscandia,var = sem.fen.precip$residuals, plottitle = "Residuals for unscaled SELM, with NDVI, Precip, Bogs")

Figure 8: Map of residuals for SELM model with NDVI, log Precipitation and Bog as predictors, on a) scaled and b) unscaled variables

s.sem.coef.precip <- readRDS("s.sem.coef.precip.rds")

scaled.plot.precip <- ggplot(s.sem.coef.precip)+geom_col(aes(x=Estimate,y=Parameter), width = 0.5, fill = "#d73027")+
  geom_errorbar(aes(y=Parameter,xmin = Estimate-std.error,xmax=Estimate+std.error), width = 0.5)+
  labs(x = "Scaled estimates",y = "Predictor")+
  theme_bw(base_size = 8)+
  theme(legend.position = "bottom")
sem.effectsize.precip <- readRDS("sem.effectsize.precip.rds")

effectsize.plot.precip <- ggplot(sem.effectsize.precip)+
  geom_col(aes(x=Percent,y=Parameter), width = 0.5, fill = "#4575b4")+
  geom_errorbar(aes(y=Parameter,xmin = Percent-(100*std.error),xmax=Percent+(100*std.error)), position = "dodge", width = 0.5)+
  labs(x = "Effect size in %", y = "")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
compare.precip <- cowplot::plot_grid(plotlist = list(scaled.plot.3,effectsize.plot.3),ncol=2)
save_plot(plot = compare.precip, filename = "plot.sem.fen.precip.png" )
print(compare.precip)

Figure 9: Comparison of scaled estimates and effect sizes from the limited SELM model with NDVI, log Precipitation and Bog as predictors.

5.5 Spatial error linear model on Fennoscandian datasets with NDVI, TNdep and logPrecip

k.neigh.w <- readRDS("k.neigh.w.Rdata")
fennoscandia <- readRDS("fennoscandia.rds")

s.fm <- s.logTOC~s.NDVI+s.logPrecip+s.TNdep
s.sem.fen.tndep <- errorsarlm(s.fm,fennoscandia,k.neigh.w)
saveRDS(s.sem.fen.tndep,"s.sem.fen.tndep.rds")

fm <- logTOC~NDVI+logPrecip+TNdep
sem.fen.tndep <- errorsarlm(fm,fennoscandia,k.neigh.w)
saveRDS(sem.fen.tndep,"sem.fen.tndep.rds")
s.sem.fen.tndep <- readRDS("s.sem.fen.tndep.rds")
sem.fen.tndep <- readRDS("sem.fen.tndep.rds")
k.neigh.w <- readRDS("k.neigh.w.Rdata")

s.sem.coef.tndep <- s.sem.fen.tndep$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
s.sem.coef.tndep <- s.sem.coef.tndep %>% rownames_to_column(var = "Parameter")
s.sem.coef.tndep$std.error <- s.sem.fen.tndep$rest.se[-1]
saveRDS(s.sem.coef.tndep, "s.sem.coef.tndep.rds")

sem.coef.tndep <- sem.fen.tndep$coefficients[-1] %>% as.data.frame() %>% 
              setNames("Estimate")
sem.coef.tndep <- sem.coef.tndep %>% rownames_to_column(var = "Parameter")
sem.coef.tndep$std.error <- sem.fen.tndep$rest.se[-1]

sem.effectsize.fen.tndep <- compute.effect.size(sem.coef.tndep)
saveRDS(sem.effectsize.fen.tndep, "sem.effectsize.fen.tndep.rds")


n <- length(sem.coef.tndep$Parameter)

s.sem.AIC.tndep <- 2*(n-1)-2*s.sem.fen.tndep$LL
s.moran.I.res.sem.tndep <- moran.test(s.sem.fen.tndep$residuals, k.neigh.w, alternative = "two.sided")

sem.AIC.tndep <- 2*(n-1)-2*sem.fen.tndep$LL
moran.I.res.sem.tndep <- moran.test(sem.fen.tndep$residuals, k.neigh.w, alternative = "two.sided")

summary.table.sem.tndep <- cbind(dplyr::select(sem.coef.tndep,Parameter),s.sem.coef.tndep[,2:3],sem.coef.tndep[,2:3],sem.effectsize.fen.tndep[,2:4]) %>%
                rbind(c("AIC",s.sem.AIC.tndep,"",sem.AIC.tndep,"","","","")) %>%
                rbind(c("Moran'I res", s.moran.I.res.sem.tndep$estimate[[1]],"",moran.I.res.sem.tndep$estimate[[1]], "","","",""))
summary.table.sem.tndep[,-1] <- apply(summary.table.sem.tndep[,-1],2,as.numeric)

saveRDS(summary.table.sem.tndep, "summary.table.sem.tndep.rds")

Table 4: Summary table for SELM model with 3 predictors: NDVI, TNdep and logPrecip

summary.table.sem.tndep <- readRDS("summary.table.sem.tndep.rds")
n <- 3

knitr::kable(summary.table.sem.tndep,digits = 3, caption = "SELM with 3 predictors (NDVI, TNdep, Precipitation") %>%
  add_header_above(c("Model"=1,"Scaled SELM" = 2,"SELM" = 2, "Effect size" = 3),italic = T) %>%
  kable_styling(bootstrap_options = "bordered") %>%
  column_spec(column = 8, bold = T) %>% 
  group_rows(group_label = "Model Evaluation", start_row = n+1,end_row=n+2)
SELM with 3 predictors (NDVI, TNdep, Precipitation
Model
Scaled SELM
SELM
Effect size
Parameter Estimate std.error Estimate std.error Estimate std.error Percent
NDVI 0.327 0.012 2.478 0.090 1.025 0.002 2.509
logPrecip -0.420 0.032 -1.166 0.089 0.988 0.001 -1.153
TNdep 0.213 0.027 0.000 0.000 1.012 0.001 1.156
Model Evaluation
AIC 7916.082 6348.003
Moran’I res 0.012 0.012 
s.sem.fen.tndep <- readRDS("s.sem.fen.tndep.rds")
sem.fen.tndep <- readRDS("sem.fen.tndep.rds")

f_plotspatial(data = fennoscandia,var = s.sem.fen.tndep$residuals, plottitle = "Residuals for scaled SEM, with NDVI, Precipitation, TNdep")

f_plotspatial(data = fennoscandia,var = sem.fen.tndep$residuals, plottitle = "Residuals for unscaled SELM, with NDVI, Precipitation, TNdep")

Figure 8: Map of residuals for SELM model with NDVI, log Precipitation and Bog as predictors, on a) scaled and b) unscaled variables

s.sem.coef.tndep <- readRDS("s.sem.coef.tndep.rds")

scaled.plot.tndep <- ggplot(s.sem.coef.tndep)+geom_col(aes(x=Estimate,y=Parameter), width = 0.5, fill = "#d73027")+
  geom_errorbar(aes(y=Parameter,xmin = Estimate-std.error,xmax=Estimate+std.error), width = 0.5)+
  labs(x = "Scaled estimates",y = "Predictor")+
  theme_bw(base_size = 8)+
  theme(legend.position = "bottom")
sem.effectsize.fen.tndep <- readRDS("sem.effectsize.fen.tndep.rds")

effectsize.plot.tndep <- ggplot(sem.effectsize.fen.tndep)+
  geom_col(aes(x=Percent,y=Parameter), width = 0.5, fill = "#4575b4")+
  geom_errorbar(aes(y=Parameter,xmin = Percent-(100*std.error),xmax=Percent+(100*std.error)), position = "dodge", width = 0.5)+
  labs(x = "Effect size in %", y = "")+
  theme_bw(base_size = 8)+theme(legend.position = "bottom")
compare.tndep <- cowplot::plot_grid(plotlist = list(scaled.plot.tndep,effectsize.plot.tndep),ncol=2)
save_plot(plot = compare.tndep, filename = "compare.fen.tndep.png" )
print(compare.tndep)

Figure 10: Comparison of scaled estimates and effect sizes from the SELM model with NDVI, TNdep and log Precipitation as predictors.

knitr::knit_exit()
CORDEX. 2021. “CORDEX — vERC.” https://portal.enes.org/data/enes-model-data/cordex.
Detsch, Florian. 2021. “Download and Process GIMMS NDVI3g Data [R package gimms version 1.2.1].” Comprehensive R Archive Network (CRAN). https://cran.r-project.org/package=gimms.
Henriksen, Arne, Brit Lisa Skjelvåle, Jaakko Mannio, Anders Wilander, Ron Harriman, Chris Curtis, Jens Peder Jensen, Erik Fjeld, and Tatyana Moiseenko. 1998. “Northern European Lake Survey 1995: Finland, Norway, Sweden, Denmark, Russian Kola, Russian Karelia, Scotland and Wales.” Royal Swedish Academy of Sciences. https://www.jstor.org/stable/4314692.
Hijmans, Robert J, Jacob Van Etten, Joe Cheng, Michael Sumner, Matteo Mattiuzzi, Jonathan A Greenberg, Oscar Perpinan, et al. 2018. “Package ’raster’ Type Package Title Geographic Data Analysis and Modeling.” https://github.com/rspatial/raster/issues/ https://cran.r-project.org/package=raster.
Pebesma, E. 2021. “Simple Features for R.” https://orcid.org/0000-0002-9415-4582 https://orcid.org/0000-0002-9415-4582%0Ahttps://orcid.org/0000-0003-2392-6140.
Tucker, Compton J, Jorge E Pinzon, Molly E Brown, Daniel A Slayback, Edwin W Pak, Robert Mahoney, Eric F Vermote, and Nazmi El Saleous. 2005. “An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data.” International Journal of Remote Sensing 26 (20): 4485–98. https://doi.org/10.1080/01431160500168686.